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Finding the ground state of a fermionic Hamiltonian using quantum Monte Carlo is a very diffi- 
cult problem, due to the Fermi sign problem. While still scaling exponentially, full configuration- 
interaction Monte Carlo (FCI-QMC) 1 mitigates some of the exponential variance by allowing anni- 
hilation of noise - whenever two walkers arrive at the same configuration with opposite signs, they 
are removed from the simulation. While FCI-QMC has been quite successful for quantum chemistry 
problems 1 ' 2 , its application to problems in condensed systems has been limited 3,4 . In this paper, we 
apply FCI-QMC to the Fermi polaron problem, which provides an ideal test-bed for improving the 
algorithm. In its simplest form, FCI-QMC is unstable for even a fairly small system sizes. However, 
with a series of algorithmic improvements, we are able to significantly increase its effectiveness. We 
modify fixed node QMC to work in these systems, introduce a well chosen importance sampled trial 
wave function, a partial node approximation, and a variant of released node. Finally, we develop a 
way to perform FCI-QMC directly in the thermodynamic limit 

PACS numbers: 05.10.Ln 

Full configuration-interaction quantum Monte Carlo (FCI-QMC) 1 is a method for finding the ground state of a 
fermionic Hamiltonian H . Starting from some state with non-zero ground state overlap, FCI-QMC stochastically 
performs imaginary time propagation \tpp) = e~P H \ipi). Working in a second quantized formalism, where anti- 
symmetry enters in the off-diagonal components of H, FCI-QMC formally yields the fermionic ground state \ipo) for 
sufficiently large (3. 

The wavefunction at any instant in imaginary time is sampled as a set of N w walkers, where each walker is a signed 
element of some many-body basis. Breaking up e~^ H into N small time steps, e~ im = e -^/ iv e -W iv . . . e -P H / N ~ 
(1 — tH) n , where r = /3/N, one applies the operator Ui = 1 — tH stochastically 5 . 

The fermion sign problem manifests in the negative off-diagonal elements of U\ , whose sign is set by Fermi statistics 6 . 
In general, there will be loops such that a walker starts in some basis state \D), hops around the loop via the off- 
diagonal elements of U±, and then returns to \D) with the opposite sign. If these sign-violating loops exist for a given 
Hamiltonian and basis, they are said to have a sign problem. 

FCI-QMC attempts to mitigate the sign problem by annihilation; any time two walkers end up on \D) with opposite 
signs, they are both removed from the simulation. It relies on efficient annihilation to prevent walkers with the incorrect 
sign from propagating as noise. In principle one can always pick a very large N w to achieve efficient annihilation, but 
the necessary number of walkers increases significantly with basis size. The efficiency of FCI-QMC has been explored 
in a few systems, including all-electron molecules 1 , the homogeneous electron gas 3 , and the Hubbard model 4 . 

In this paper we develop extensions to FCI-QMC and apply these new developments to the Fermi polaron 7 . The 
polaron is a classic problem in condensed matter physics, as it is one of the simplest problems that shows strongly- 
interacting many body effects. We will be particularly interested in the Fermi polaron as seen in cold atomic gases 8 , 
whose ground state properties are well understood through a combination of analytical 7 and numerical 9 methods. As 
such, the polaron problem serves as a nice test-bed for testing and improving the FCI-QMC algorithm in a condensed 
matter setting. 

In the remainder of this paper, we introduce the polaron problem and discuss how we calculate its properties by 
improving on the FCI-QMC algorithm. In particular, we develop approaches which start with a sign-free Hamiltonian 
and reintroduce these signs so as to achieve an exact answer. In addition, we develop a way to do FCI-QMC in the 
thermodynamic limit. 

In section I, we introduce the polaron. Then, for a test value of 10 6 walkers, we find that the naive implementation 
of FCI-QMC works for a basis restricted to only a single excitation on the non-interacting ground state, but fails for 
cases with additional excitations. We proceed to describe a trial wavefunction (section III) inspired by variational 
solutions to the polaron problem 10 and introduce a variant of FCI-QMC which performs importance sampling. This 
increases the accessible parameter regime, but not sufficiently. Next, we introduce a controlled approximation (section 
IV) to attenuate the sign problem, which we call the partial node approximation. We extrapolate observables deduced 
from partial node simulations to the exact answers. We then combine release node techniques (section V) with the 
partial node approximation as an alternative way to find the true ground state observables. The partial node energies, 
although extrapolable to the correct result, are not individually variational. At the cost of a time-step error, we show 
(section VI) how to recover the variational guarantee by performing a more typical fixed-node approach. To accomplish 
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this, we overcome the non-trivial obstacle that a single basis element is connected by the Hamiltonian to an extremely 
large number of other basis elements. We discover that all three approaches for computing the ground state properties 
are consistent with each other and in agreement with other analytical and numerical work (section VII). Finally, we 
note that extrapolation to the thermodynamic limit is made difficult by the presence of shell effects. To avoid shell 
effects, we introduce an extension of the FCI-QMC algorithm that allows for working directly in the thermodynamic 
limit (section VIII). We compare and contrast our thermodynamic limit solution with that of diagrammatic Monte 
Carlo, and comment on future possibilities for utilizing our improvements to FCI-QMC in the field of condensed 
matter physics. 



I. THE POLARON PROBLEM 



The polaron problem we consider is defined for a three-dimensional system of two fermion species, denoted spin 
up and spin down. Starting with a non-interacting Fermi sea of spin up particles, a single spin down "impurity" is 
added that is able to interact with the sea of spin ups. Experimentally, near a broad s-wave Feshbach resonance the 
fermions interact by a short-range potential with scattering length a. The momentum space Hamiltonian is 7 

h = £kci a c k<J + ~ c l+ q .,A-q,i c p^ ck ^ ' w 

ka kpq 

where e k = k 2 j (2m), V is the volume of the system, and m is the mass of the particles. We work in units where H = 1. 

We regularize the short-range interaction by introducing a high-momentum cutoff A on the spin-ups. The interaction 
strength g is cutoff-dependent 11 : 

^=8^-4^ (2) 

We solve this Hamiltonian for a cube of side length L = V 1 / 3 with periodic boundary conditions. With N spin up 
particles, the Fermi wavevector and energy are kp = (6tt 2 N /V) 1 / 2 and Ep = k F /(2m) respectively. 
In the non- interacting limit (a = 0), the ground state is an undressed polaron, 

|X> ) = |FS t , 4 > , (3) 

where |FS^,0^) denotes a Fermi sea of spin-up particles with a single spin down at zero momentum. Upon turning 
on a weak attractive interaction (small negative a), the polaron remains a well-defined quasiparticle with non-zero 
quasiparticle residue Z = \(4>o\Dq)\ 2 , where IV'o) is the ground state of the interacting Hamiltonian. Note that Z can 
be written as the ground state expectation (ipQ\Vo\ipo) , where Vq = |Do)(-Do| projects onto the undressed polaron. 

The polaron problem has been most thoroughly investigated for the case where the scattering is unitary limited, 
(kpa)~ x = 0. Chevy found that the ground state for these parameters is well-described by a variational wave- 
function that includes all single particle-hole excitations on the spin up Fermi sea, giving a ground state energy of 
— 0.6066£'f 7 ' 10 - These results were later extended by Combescot and Mora 10 , who introduced a procedure for calcu- 
lating the variational ground state energy at an arbitrary number of particle-hole excitations and proceeded to solve 
for the two particle-hole pair variational energy, — 0.6156-Ef. 

The Fermi polaron has also been approached via Monte Carlo techniques, most notably diagrammatic Monte Carlo 9 . 
Monte Carlo results for the ground state energy agree well with the variational solutions. Given the agreement of 
the ground state polaron energy between variational expansions, Monte Carlo, and experiment 8 , the energetics of the 
polaron at unitarity are a good test case for improvements to the FCI-QMC algorithm. 

Off of unitary on the BEC side (a" 1 > 0), the polaron eventually becomes unstable to formation of a tightly-bound 
molecule. This is theoretically predicted to occur as a first order "phase transition" 11 near (fcpa) -1 = 0.9. However, 
experimentally the quasiparticle residue vanishes at (fcjfo)" 1 « 0.75 8 . It is believed that this difference is due to the 
small but finite density of the spin down particles in the experiment; we will address this issue in sec. VII. 



II. FCI-QMC APPLIED TO THE POLARON 



Given the effectiveness of the variational expansions in describing the polaron, we choose as our basis all momentum 
space determinants with at most M particle-hole pairs dressing the spin-up Fermi sea. The Chevy ansatz 7 is then, 
for example, the exact solution for M = 1. In this basis, all off-diagonal elements of H come from the interaction and 
are given by the constant value g/V, with their sign set by Fermi statistics. We fix the number of spin-up particles N 
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TABLE I: Basis size and variational parameters a and /3 used for different Hamiltonians defined by [A, M, N, l/(fe/a)], as 
described in the text. 





FIG. 1: (color online) (a) Polaron energy traces at (kpa) -1 = 0.9 for importance sampled FCI-QMC with N = 33 spin up 
particles, momentum cutoff K/kp = 20, M = 1 (dashed blue), and M = 2 (solid red). Discontinuities in the M — 2 data are 
regions where the denominator of the energy metric (5) passes through zero, (b) Histogram showing the momentum distribution 
of the particle excitations for the M = 1 data in part (a). Momenta are uniformly occupied from the edge of the Fermi sea 
(k = kp) up to the cutoff (k = 20kp). (c) Graphical representation of sign issues for the polaron. One can show that single 
particle-hole excitations all have positive signs in the exact ground state, but determinants with two particle-hole pairs are 
connected to them by signs that are half positive and half negative. 



and momentum cutoff A. Note that we eventually want to take the physical limit N — > oo and A — > oo, for which the 
basis size becomes infinite. We work in the sector with zero total momentum, which is always valid for a polaron-likc 
quasiparticle 11 . N will be limited to 7, 19, 27, ... to ensure that the spin up Fermi sea is a closed shell. 

The walkers are initialized in the undressed polaron state with positive sign, = \Dq). Each walker is labeled by 
the pair {S w , \D W )}, corresponding to its sign and determinant respectively. We will write general determinants as 

\D) = |ki,k 2 ,...,k M ;qi,...,qM) = 4 ± •••cj CM c qi •••c qM |D ) , (4) 

where k's label the particle excitations, q's label the holes, kx < k 2 < . . . < kjyi, and qi < . . . < qjyi for an arbitrary 
but fixed momentum ordering. The values of ki and qi are restricted to momenta on a discrete grid set by the size of 
the box. The energy metric in this basis is 

(D \H\j> w (P)) = T, w S w (Do\H\D w ) 
[P) (A>|Vv,(/?)> Y, W S W (D \D W ) ' {0) 

where the zero of energy is defined as (Dq\H\Dq) = 0. For sufficiently large /3, the average of E(j3) converges to the 
fermionic ground state energy. 

At unitarity with N = 33 and A = 20, FCI-QMC is able to find the ground state energy for M = 1, where the 
sign problem is weak (fig. la). For M = 2 however, the sign problem prevents convergence with fixed number of 
walkers N w = 10 6 (fig. la). The M = 2 sign problem can be understood by starting from a determinant with 
two particle-hole pairs, I-D2) = |ki, k 2 ; qi, <ja) • Off-diagonal matrix elements allow hopping from \D 2 ) to four single 
particle-hole determinants (see fig. lc). Of these, two come with positive sign, two with negative sign, and all have 
the same weight g/V. So the determinants with two particle-hole pairs rain down nearly-random signs on the single 
particle-hole shell, which one can show should all be positive in the ground state. This results in a large sign problem 
unless the random signs are efficiently annihilated. 
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We emphasize that our approach to the Fermi polaron differs from previous problems where FCI-QMC has been 
applied because 1) we cut off the Hamiltonian after a fixed number of particle-hole pairs, 2) we have a single particle 
basis that is extremely large , and 3) almost all determinants are important. The last fact arises because the interaction 
is short-range in real space, i.e. long-range in momentum space. Fig. lb shows a histogram of the walkers in the first 
particle-hole shell for the M = 1 ground state. The particle excitations are evenly occupied from the Fermi sea all the 
way up to the cutoff, in contrast to situations 1 where the ground state consists of a few large-weight determinants; 
this makes annihilation more difficult. 

FCI-QMC for quantum chemistry problems minimizes these issues through an intelligent choice of single particle 
orbitals. However, for the Fermi polaron, there is no clear choice of basis that minimizes the impact of unimportant 
determinants. In the next section we show that using importance sampling can enhance the quality of the algorithm. 
This not only improves the variance, as usual in importance sampling, but also improves the annihilation properties 
by favoring certain determinants. 



III. IMPORTANCE SAMPLING AND TRIAL WAVE-FUNCTION 



Importance sampling is by now a standard Monte Carlo technique that has been used to decrease statistical noise. 
Importance sampling consists of sampling walker \D) from a probability distribution proportional to \(D\iPt) (D\^jq)\ 
instead of K-DlV'o)!; where \iPt) is a trial wavefunction that we choose. This is implemented by re-weighting the 
off-diagonal moves, and can be thought of as simply acting with the (non-Hermitian) effective "Hamiltonian" H{ s 
given by 



The energy metric becomes 



(D>\H is \D) = (D>\H\D) { ^^± . (6) 



F(B] = E w S w (Do\H\D w )/(jj T \D w } 
[> J2 W S W (D \D W )/^ T \D W ) ■ {) 



We also calculate quasiparticle residue through the method of mixed estimators 12 , 

Z « 2(ifr\V \ifo) - (V>t|Wt) ■ (8) 

As His is not Hermitian, we must be careful to only act with it on kets. 

Choosing an appropriate trial wavefunction plays a major role in the success of importance sampling. A naive guess, 
the free fermion wavefunction, would give |V>t) = |A)); this would not allow sampling of any particle-hole excitations. 
Therefore, we must construct some IV't), which in general we would like to be as close as possible to the ground state 

To motivate our choice of \i/>t)> consider the Schrodinger equation H\tpo) — Eo\ipo)- H H = T + V, where the kinetic 
term T is diagonal in the momentum basis, then 

(J| ^ = - (D|f|D)- £o g (J|m ' >(iyW - <9) 

For \D) with n particle-hole pairs, the interaction term connects it to determinants \D') with n — 1, n, and n + 1 
particle-hole pairs. The strength of the interaction is g/V with sign Sd> = sgn((D\V\D'}). While these coupled linear 
equations are generally hard to solve, we can make the simplifying approximation of only considering \D') with n — 1 
particle-hole pairs in constructing \iI>t)- To allow some freedom, we will define variational coefficients a and j3 to 
replace g and Eq respectively, giving 

(DM = - J- S^jD'^r) , (10) 

(D n \T\D n )-a D' n _ 1 s.t. 

where \D n ) denotes a determinant with n particle-hole pairs. This gives a recursive definition for \i/jt) which, since 
normalization is irrelevant, we seed from (Dq\iPt) — 1- The parameter a is chosen to be approximately the energy. 
Then, for each shell of n particle-hole pairs, we calculate Z n — (ipr\'Pn\'<l'T) using variational Monte Carlo, where V n 
projects onto the subspace of determinants with n excitations. j3 is chosen to minimize the maximum deviation of Z n 
from its average value, 1/(M + 1). 
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FIG. 2: (color online) Polaron energy (a) and quasiparticle residue (b) for N = 33, (Aifo) -1 = 0, A/fcp = 20, with M = 1 
to 4 particle- hole excitations (to improve statistics we use 7 = 0, so energies are a variational upper bound). Panel (c) shows 
energies upon breaking off of unitarity. Error bars are all computed for runs of fixed imaginary time A/3 = 35/ Ef and therefore 
can be reasonably compared with each other. No data is shown for (fc^a) 1 > 0.15 with M — 2 because the simulation failed 
to converge a sign structure in that regime. 



There is one minor issue with this ipr- due to the discreteness of the Hilbert space, eq. (10) will potentially yield 
(D\iPt) — for some small set of determinants \D), which become increasingly scarce as N — > 00. This prevents any 
weight from being put on these determinants, introducing a bias. To avoid this bias, we introduce a third parameter 
7 such that if eq. (10) yields |(.D|-0t)| < 7, we set (D\ipx) 7- The biased simulation 7 = gives a variational upper 
bound on the energy, which should be fairly close to correct. As an example, at unitarity for the relatively small basis 
with N = 7, A = lOJfejr, and M = 3, the difference in energy is -0.6279(3)^ for 7 = versus -0.6310(17).E F for 
7 = 10 -4 . We use 7 = 10~ 4 throughout the remainder of this paper, unless otherwise specified. 

Importance sampling greatly improves the effectiveness of FCI-QMC. For example, at unitarity FCI-QMC with 
importance sampling is able to solve the polaron ground state for A = 20kp, N — 33, and M = 4 - corresponding 
to a basis size of 2.03 x 10 21 - using only 10 6 walkers. Solutions for the polaron energy and quasiparticle residue at 
unitarity are shown in fig. 2a,b. 

We next push to positive values of a -1 , where for (kpa) -1 > 0.9 the polaron is expected to become unstable 
to formation of a molecule. As the interaction strength X/(kpa) is increased, the polaron becomes more strongly 
correlated. The weight of the wavefunction is pushed to higher momenta, making annihilation more difficult. As seen 
in fig. 2c, for M = 2 particle-hole pairs, the error bar gradually increases as we push off of unitarity. Eventually, for 
l/(kpa) > 0.15, FCI-QMC with importance sampling fails to find the ground state. To proceed, we now introduce 
modifications that allow us to further stabilize the algorithm. 



IV. PARTIAL NODE APPROXIMATION 



Fixed node quantum Monte Carlo is a method that removes the exponential cost of solving a fermionic Hamiltonian 
by approximating the Hamiltonian with a related one that has no sign problem. In this section we instead find a way 
to extrapolate to the exact solution of our (restricted basis) Hamiltonian using a fixed node- inspired starting point. 
In fixed node algorithms, such as the lattice formulation of van Bemmel et. al. 13 , one must specify a "correct" sign 
for each determinant. For simplicity we choose to use \iPt) as determining our sign structure. Given this choice, an 
off-diagonal element (D'\H- iS \D) will be called sign-violating if (D'\Hi a \D) > 0, since in this case applying —rH m to 
hop from D to D' would flip the sign, violating the sign structure of \ipr)- 

The fixed node Hamiltonian Hf n is given by 

(D'\H ln \D) = I ° if Sign vi ° L ^ , where D ? D' 

V 1 ' ' \(D'\H is \D) if not s.v. (n.s.v.) ^ 

(D\H ln \D) = (D\H is \D)+ ]T (D'\H is \D) . (11) 

D'S.V. 

Note that the sign-violating off-diagonal matrix elements have been removed and "dumped" onto the diagonal. This 
Hf n is guaranteed to give a variational upper bound on the ground state energy, yielding the correct ground state 
energy if \ip T ) =\ip ). 

Lattice fixed node was originally designed for real space lattices where there are very few off-diagonal elements 
connected to any given configuration. For our Hamiltonian, there will in general be many sign-violating matrix 
elements connected to a given determinant \D) (of order > 10 7 for our parameters), and the sums required to modify 
the diagonals of Hf n become analytically intractable. We therefore initially work with a modified Hamiltonian H' fn - 
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FIG. 3: (color online) (a) Partial node FCI-QMC energies, both with (red triangles) and without (blue dots) diagonal dumping 
(see text). Note that the e pn = result, which is computationally feasible for M = 1, is just the solution for the exact 
Hamiltonian H using importance sampling. A linear extrapolation of the partial node energies (green star) agrees well with the 
e pn = solution, (b) Partial node energies for M = 2, where the sign problem prevents a solution of the exact Hamiltonian. 
A quadratic extrapolation agrees well with the fully-signed solution (black square), which is obtained via release node QMC 
initialized from the ground state walkers of partial node FCI-QMC with e pn = 0.4. (c) Average energy for 100 release node 
traces used in panel (b). Arrows indicate the region used in calculating statistics. Error bars show standard error across trials. 
The inset shows a similar release node trace in which the initial condition is intentionally far from correct, to make the decay in 
energy more apparent, (d) Flip rate of the denominator of the energy metric as a function of e pn (blue dots). For comparison, 
the red triangle shows the flip rate without importance sampling and the black dashed line shows the inverse of the average 
the Monte Carlo correlation time for our simulation. All data are for l/(kpa) = 0.9, A = 20/cf, and N = 33. (a)-(c) use 10 6 
walkers, while (d) uses 10 4 . 



which we refer as the fixed node Hamiltonian without diagonal dumping - given by 

(D>\HL\D) = I if Sign vioL , where D ? D' 

V 1 fnl 1 \ (D'\H is \D) if not s.v. ' r 

(D\H' ln \D) = (D\H is \D) , (12) 

whose ground state energy is no longer guaranteed to be a variational upper bound on that of H . In section VI we 
will describe the steps needed to reintroduce the upper bound. 

Due to the presence of annihilation in FCI-QMC, we are able to introduce a less-severe approximation, which we 
call the partial node approximation. In partial node, we interpolate between the exact and fixed-node Hamiltonians 

-ffpn = e pn-fffn + (1 — £pn)-ffis • (13) 

Surprisingly, we find that this algorithm works up to a very high fraction of the original sign violating terms (small 
e pn ). To obtain the exact answer, we extrapolate to the original Hamiltonian (e pn ~ 0), heuristically using a quadratic 
fit (see fig. 3). 

To quantify the improvement gained by using the partial node approximation, we measure the average flip rate of 
the sign of the walkers on the reference determinant Dq (which typically should contain approximately 30% of the 
walkers for a given snapshot) . The flipping of this sign is a signature of an instability of the algorithm coming from the 
sign problem, and the average flip rate must be significantly below the inverse correlation time for useful data to be 
garnered from the simulation. Fig. 3d shows the average flip rate using the partial node approximation. Notice that 
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£pn = is close to the inverse correlation time, presenting a problem, but by e pn = 0.4 the flip rate has significantly 
decreased. 



V. RELEASE NODE 

In FCI-QMC without annihilation, the sign problem manifests as a variance that scales exponentially with both 
P and the energy difference between the fermionic and bosonic ground states 6 . This exponential scaling in j3 holds 
even in the presence of annihilation at finite walker number, as annihilation can be thought of as simply reducing the 
fermion-boson ground state gap. Therefore, minimizing the total /3 needed to reach the ground state wave-function 
with a certain fidelity also minimizes the effect of the sign problem. Thus, one can improve the convergence time of 
the importance-sampled FCI-QMC algorithm (e pn = 0) by starting very close to the ground state. 

One approach to starting near the true ground state is to use the ground state of H pn with a small value of e pn . We 
aren't able to explicitly represent this wavefunction, but we can prepare it stochastically via partial node FCI-QMC. 
The set of N w walkers thus prepared is then allowed to evolve under the true Hamiltonian for fixed imaginary time 
/3. When the wavefunction is "released" to evolve under the exact Hamiltonian, it quickly relaxes to the exact ground 
state; this is what we call release-node FCI-QMC 14 . 

In release node, we use the same energy metric as in importance sampled FCI-QMC. In this case, since we have 
now propagated with two non-commuting Hamiltonians (H pn and Hi S ), one can show that the energy metric only 
becomes meaningful when the wavefunction reaches the ground state, at which point it gives the ground state energy. 
As seen in figure 3, the energy converges to the ground state energy before the errors blow up if we start close enough 
to the ground state. The release node energies agree well with the partial node extrapolation (fig. 3), which provides 
a check of the validity of our extrapolation. 

VI. DIAGONAL DUMPING 

The partial node approximation that we have discussed thus far does not give a variational upper bound on the 
energy, since we removed the diagonal dumping to get the fixed node Hamiltonian H' {n (12). To restore the variational 
upper bound, we now discuss how one can put the effects of the diagonal dumping back into the partial node algorithm. 

While the sum involved in diagonal dumping is in general not analytically tractable, it can be done stochastically 
in much the same way as off-diagonal spawning. During the diagonal create/kill step 1 we would like to stochastically 
apply 

t/dia g = 1 - T(D\H is \D) - r ^ (D'\H is \D) (14) 

D'S.V. 

to a walker on determinant \D). The first two terms can be done exactly, while the sum - denoted AK - can be 
sampled as follows. 

1. From \D) pick another determinant \D') according to some normalized probability function p gcri {D'\D). 

2. If (D'\H is \D) is not sign- violating, then AK = 0. 

3. Otherwise, AK = (D'\H is \D) /p gcn (D'\D). 

However, as AK is now potentially very large, multiplying \D) by weight t/diag has the potential to be disastrous if 
Kiiag <C — 1. Therefore, we find it necessary when stochastically dumping the diagonal to also use the approximation 

1 - r(D\H is \D) - tAK » e -r((D\H is \D)-AK) (15) 

This introduces a time step error which one must extrapolate to zero; empirically a quadratic extrapolation for 
t < 10~ 3 appears to be quite effective. Partial node results with stochastic dumping of the diagonal are shown in fig. 
3. 



VII. RESULTS 

Using partial node FCI-QMC followed by release node, we are able to determine values for the energy and quasi- 
particle residue of the Fermi-polaron on the BEC side of the interaction, which are shown in fig. 4 for A = 20kp and 
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FIG. 4: (color online) Ground state energy (a) and quasiparticle residue (b) for the polaron off of unitarity at A = 20kp 
and N = 33. Energies are compared against the TV — > oo, A — > oo results from diagrammatic MC 9 , while Z is compared to 
experimental measurements 8 and the analytical variational result at M = l 7 . Deviations of Z from the analytic results likely 
come from a combination of finite size effects and error in the mixed estimator, (c) and (d) show finite size effects seen in 
attempting to extrapolate the energy in 1/A, with (feFfl) -1 = 0, = 33, M = 1(c) and 2 (c,d). In (d), we compare energy 
against single particle basis size, i.e. the number of available spin up momenta < k < A. 

N = 33. The energies are compared to polaron and molecule energies from diagrammatic Monte Carlo 9 . Our results 
match well for small values of l/(kpa), while for larger values of l/(/cpa) our energies deviate from the diagrammatic 
results. This is a consequence of working at fixed A = 20kp, as opposed to diagrammatic Monte Carlo, which is done 
in the limit A — > oo. Our results match with the diagrammatic solution after extrapolating to A = oo. 

Figure 4b shows the quasiparticle residue Z compared to experimental results 8 , where there is a small but finite 
density of spin-down atoms. We find that, as with the energy, an M — 1 variational expansion provides a good estimate 
of Z. Therefore, as noted elsewhere 8 , our theoretical model differs significantly from the experiment, possibly as a 
result of the finite spin-down density used experimentally. 

In extrapolating our results to the physical limit, A — > oo, we encountered an unexpected problem. Using a linear 
fit over a wide range of 1/A (fig. 4a), we find that individual data points have error well outside the line, with no 
discernible pattern. Zooming into a very small region of A we discovered the reason behind this: small fractional 
increases in basis size as new shells become available cause large jumps in the ground state energy. A similar issue 
occurs when attempting to extrapolate in particle number N. We nevertheless attempt such an extrapolation in fig. 
4 a and b. As a result, the size of the error bars reflects not the accuracy of the data points at individual values of 
A, but rather these inherent shell effects. 

Shell effects will vanish in the limit N — > oo. Therefore, we next address the possibility of extending FCI-QMC 
into the thermodynamic limit. 

VIII. FCI-QMC IN THE THERMODYNAMIC LIMIT 

All the quantum Monte Carlo simulations described so far have been done for a finite number of particles. In this 
section, we describe how to modify FCI-QMC to work directly in the thermodynamic limit (TDL, N —> oo). In the 
thermodynamic limit, the accessible momenta span a continuous set of k-points instead of being limited to a discrete 
grid. Additionally, it is important to work in a representation where instead of enumerating the momenta for all 
N = oo particles, we instead store only their excitations above a known state, in this case \Dq). Crucially, in order 
for the spectrum of the Hamiltonian to remain bounded, allowing application of 1 — tH instead of eT rli , the TDL 
only works with a cutoff M on the number of excitations allowed. We note that there are modifications to FCI-QMC 
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FIG. 5: (color online) Ground state polaron energy in the TDL (a) with M = 1, (fc.Fa) _1 = as a function of A (blue dots). 
A linear extrapolation in 1/A (green star) agrees with the exact ground state energy with M — 1, calculated by Chevy's 
variational ansatz . (b) Fixed node (e pn = 1) energies for M = 2 compared against the variational ansatz . 

that can allow it to work in continuous time, i.e. allow for application of e~ rH , which we discuss in appendix A. This 
requirement is similar in spirit to the diagrammatic Monte Carlo method of imposing a cutoff on diagram order, in 
an attempt to avoid a divergence. 

As momentum space is no longer discretized in the TDL, one might assume that annihilation is no longer possible, 
and therefore that the algorithm is doomed to fail. However, there is one key exception: due to the discrete choice of 
\Dq), annihilation can still occur at that one determinant. 

Given this, the TDL algorithm is in practice nearly identical to the algorithm with finite N - in fact, one can think 
of the TDL algorithm as just the limit of the finite algorithm for larger and larger N. We have carefully, though 
implicitly, defined \iPt) such that all factors of V cancel out when determining relevant quantities, such as spawning 
probability or energy. This is important because taking the limit N — > oo while at fixed spin-up density (constant 
Uf) means taking V — > oo as well. 

The TDL algorithm is remarkably effective at finding the ground state with M — 1 at unitarity, where the sign 
problem is weak. These energies at various values of A are shown in fig. 5a, where an extrapolation to A = oo is 
much smoother than for finite N due to the absence of shell effects. 

For the sign-problem-heavy M = 2 case, we instead show data in the sign-free fixed node limit (e pn = 1), using the 
diagonal dumping method described in sec. VI. This gives a good energy for the M = 2 polaron, which is a variational 
upper bound on the known energy of — 0.6156EV within error bars. We believe this latter approach of combining 
fixed node with working directly in the thermodynamic limit will have wide applicability even for the approximate 
calculations that currently dominate the fermion QMC literature. 

IX. DISCUSSION 

We have shown that the FCI-QMC algorithm can solve the Fermi polaron problem, after improving the algorithm 
with a smart importance sampled wave function, the introduction of partial and fixed node approximations, and the 
utilization of release node methods. 

We believe that our work demonstrates two main points that should be useful as FCI-QMC and its variants are 
applied to future problems in condensed matter physics. First, for strongly correlated condensed systems where many 
determinants are occupied in any single particle basis, we have demonstrated that physical insight - in the form of 
a good choice of trial wavefunction - can significantly improve the behavior of the FCI-QMC algorithm. Combining 
the improved statistics of importance sampling with sign-attenuating approximations such as partial node, we showed 
a significant increase in the effectiveness of FCI-QMC in solving the polaron problem. Finally, we showed that these 
methods can be made exact, up to statistical noise, via extrapolation or release node QMC. 

Second, we have shown that under certain conditions the FCI-QMC algorithm can be extended to work with systems 
in the thermodynamic limit. Furthermore, we have introduced a fixed node algorithm in this limit, which has histor- 
ically been an important tool for solving fermionic systems with QMC. We anticipate that these new developments 
will open a variety of problems in condensed matter physics to be approached using these new developments. 
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Appendix A: Continuous-time algorithm for FCI-QMC 

In this appendix, wc introduce a continuous-time FCI-QMC algorithm that has no fixed time step. We note that, 
while we have not yet implemented this algorithm, the number of off-diagonal terms that must be sampled to remove 
the penalty method errors (discussed below) is likely to make the algorithm significantly slower than the finite time 
step algorithms described earlier. 

Our continuous-time formulation of FCI-QMC is loosely based on a continuous-time lattice algorithm found 
elsewhere 15 . We will begin by introducing the algorithm for an arbitrary Hamiltonian H, which in general can 
be some non-Hermitian effective Hamiltonian, as found in importance sampling. We will start by assuming that all 
off-diagonal sums can be computed analytically, and later generalize this to the case where certain sums must be done 
stochastically. We generalize the existing algorithm to allow for situations where H has a sign problem, so off-diagonal 
elements of H will not be required to be negative. 

We would like to apply the propagator Up A — e~^ AH , where /3 a is now some fixed imaginary time. We refer to /3 a 
as the annihilation time, and think of applying Up A to each walker before performing annihilation. Breaking this up 
into small time intervals r, it becomes Up A = (1 — tH)(1 — tH) ■■■(! — tH). To perform continuous time QMC, we 
would like to take the r — ¥ limit of this expression. 

Consider applying U T = 1 — tH stochastically to some determinant \D) in the limit r — > 0. We define on-diagonal 
component Kl(D) and off-diagonal sum Vl(D), where 

K L (D) = (D\H\D) (Al) 
Vl(D) = \( D '\ H \ D )\ ■ 

D'^D 

We can break U\ up as 

1-tK -tV 
Ui = - (l-p 5 )H Ps = #i(l -Ps) +Vxp s , where 

1 - Ps Ps 
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1 - tK -tV 
K x = and Vi = (A2) 

1 - Ps Ps 

Then wc stochastically apply K\ with probability 1 — p s or V\ with probability p s . 

Furthermore, we would like to choose p s such that V\ simply corresponds to deterministically moving to \D') with 
probability proportional to \ {D'\H\D)\. Thus, 

]T \(D'\V 1 \D)\ = l=^ Ps = rV L (D). (A3) 

D'^D 

So the algorithm proceeds as follows: start from some determinant \D), and either apply K\ with probability 1 — p s 
or Vi with probability p s . V\ corresponds to hopping to a new determinant. K x simply multiplies \D) by a weight 

Wi = 1-tKUD) = 1-tK l (D) 
l-Ps l-rV L (D) 

V\ is chosen for the first time at step N with probability P{N\D) = (1 — p s (D)) N . If A = /?s/r, time f3s prior to 
spawning will come from the probability distribution 

P a (0s) - (1 - rV L f s/T T -^>° e -f )sV ^ D \ (A5) 

Given a choice of Ps from this distribution, the walker will also pick up a total weight W during the A— 1 non-spawning 
steps, where 

W(fis) = ttf s/t " 1 ^ e -0s[K L (D)-V L (D)] (A6) 

Note that, for a non-sign-violating Hamiltonian, the term in the exponent of W(J3s) is just the local energy. 

Therefore, for each walker we stochastically propagate for a time (3a, during which it will both hop and pick up 
weight; this can be done to all the walkers in parallel. Finally, we take all the walkers, perform annihilation, measure 
observables, and repeat. For Hamiltonians with relatively few off-diagonal terms, i.e. the real-space Hubbard model, 
this is the complete algorithm. For much larger off-diagonal sums, things become more complicated. 

The remainder of this discussion will describe our algorithm for Hamiltonians with large or infinite number of 
off-diagonal terms in the sums. There are two sums (integrals in the TDL) that we now perform stochastically: the 
diagonal dumping 

AK L (D)=e pn J2 \{D'\H is \D)\ (A7) 

D's.v. 

and the local potential energy sum 

V L (D)= J2 \(D'\H is \D)\ + (l-e pn )J2\(D'\H is \D}\. (A8) 

D'n.s.v. D's.v. 

The sums themselves are fairly straightforward. Using the same method as spawning in FCI-QMC, we have a 
method for generating determinant \D'} connected by an off-diagonal component of H to the starting determinant 
\D). If the probability to generate \D') is p gen (D'\D), then AK L for example is just the expectation of the observable 

_ 7} DD .(D'\H\D) / if not sign viol. 

O AKl = e P n— — , . . where n DD , = i b (A9) 

Pgen(-L> '\D) M if sign viol. 

A similar observable Oy L can be defined for computing Vl(D). Assume that we have sampled a total of N sum 
determinants D' to simultaneously calculate AKl and Vl with some mean O and standard error <j(0) for each 
observable. 

Formally, we could perform this sampling in the limit A sum — > oo, determine AKl and Vl exactly, and simply 
use them in the earlier procedure. However, for finite A^ sum , we want to utilize a variant of the penalty method to 
minimize the bias due to statistical uncertainty. Consider first the simpler case of applying the weight W = e -Ps&K L ^ 
Assume that the actual value of AKl is drawn from a Gaussian distribution with mean /i and width a, which are 
given by the expectation value and standard error of Oak l ■ Then the weight we should apply is simply 

/oo 
e-^ AK -G^(AK L )d(AK L ) = e -^-f^/2) (A1Q) 
-OO 
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A more complicated question is how to sample fis in an unbiased way. Again, assume that Vl is drawn from a 
Gaussian G^^{Vl)- Furthermore, assume that we have chosen N sum large enough that we don't have to worry about 
the negative Vl tail of the Gaussian, i.e. fi^> a. We want to sample from the full distribution 

POO 

p(Ps)= / p((3s\V L )G^(V L )dVL (All) 
Jo 

We can expand to second order in Vl — (J-, ignoring the first order term (whose integral vanishes). Then 

,-Psn r°° 



g 

p(Ps) oc 



l-(^) 2 (l + ^ + /3V/2) + 0(f^ Xr 
fi J \\ mu 



2 

a 



G^{V L )dV L (A12) 
(A13) 



/'- 

Normalizing and integrating this result, we find that the normalized cumulative distribution function is 
_ l . (sinhx - coshx)(2 + s 2 (x 2 + 4x + 6) 

P(x = uPs, s = a n) = 1 + ^ '\ '- A14 

os z + 2 

We can then sample from this distribution by picking z G [0, 1] at random, then finding the value of fts where 
P(x, s) = z. 

In summary, here is our algorithm for continuous partial node FCI-QMC: 

1. Start with N w walkers, each with weight W w — 1, sign S w , and determinant \D W ). At the first step, all walkers 
are initialized in \D$) with positive sign. At this point, there should be no walkers with the same determinant 
but opposite sign. 

2. Propagate each walker independently by e~^ A ^ H ~ s ^ as follows. At the start of this portion, define a variable j3 w 
for each walker, and initialize f3 w = 0. 



(a) For a walker in determinant \D), sample off-diagonal elements \D') for N sum steps to get Vl{D), <t\Vl{D)], 
AK L {D), and a[AK L (£>)]• 



(b) Sample /3s from P(x, s) as described above, where x — /3s/Vl(D) and s = (t[Vl(£>)]/Vl(Z)). 

(c) If (3 W + (3s < Pa, then move the walker to a new determinant \D'). To sample with weight proportional to 
\(D'\H\D)\, run a Metropolis algorithm for iV me t steps. The sign of \D) is (not) flipped when the matrix 
element (D'\H\D) is positive (negative). 

(d) Multiply the weight by e' 138 ^ ^ ^^, with the stochastic portion of e-^^l^l 15 } re-weighted via the 
penalty method as described above. If j3 w + fis > Pa , then use /3a — flw in place of 0s in these formulas. 

(e) Increment /3 W by /3s- If Pw < Pa, repeat from part (a). 

3. Annihilate the walkers, keeping track of their weights. For example, if five walkers are all on determinant D 
with weights W\ to W$ and signs Si to S5, let W = J2i WiSi. Then the five walkers are replaced by a single 
walker with weight \W\ and sign sgn(VF). 

4. Resample the walkers. Replace a single walker with weight W by \ W\ walkers of weight 1. Add another walker 
of weight 1 with probability W — \ W\ . 

5. Measure observables. 

6. Adjust S using standard feedback protocols 1 as desired. Then repeat from step 1. 



